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Abstract 

For the investigation of higher order Feynman integrals, potentially with tensor structure, it is 
highly desirable to have numerical methods and automated tools for dedicated, but sufficiently 
'simple' numerical approaches. We elaborate two algorithms for this purpose which may be 
applied in the Euclidean kinematical region and in d = 4 — 2e dimensions. One method uses 
Mellin-Barnes representations for the Feynman parameter representation of multi-loop Feynman 
integrals with arbitrary tensor rank. Our Mathematica package AMBRE has been extended for that 
purpose, and together with the packages MB (M. Czakon) or MBresolve (A. V. Smirnov and V. 
A. Smirnov) one may perform automatically a numerical evaluation of planar tensor Feynman 
integrals. Alternatively, one may apply sector decomposition to planar and non-planar multi-loop 
e-expanded Feynman integrals with arbitrary tensor rank. We automatized the preparations of 
Feynman integrals for an immediate application of the package sector_decomposition (C. 
Bogner and S. Weinzierl) so that one has to give only a proper definition of propagators and 
numerators. The efficiency of the two implementations, based on Mellin-Barnes representations 
and sector decompositions, is compared. The computational packages are publicly available. 
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1 Introduction 



One goal of present calculations in particle physics is reaching higher and higher precision in 
perturbation theory. Since Feynman's time we have rules allowing to build automatically the 
necessary mathematical objects. For a long term these were just the Feynman diagrams, but 
recently other approaches like the unitarity based perturbative approach get rising attention. The 
problem remains how to evaluate the complicated integrals originating in the perturbative picture 
of elementary particle interactions. Physics predicts these integrals to be defined in Minkowskian 
space-time. The integrals are multi-dimensional complex functions with a complicated singularity 
structure. Further, they have to be regularized by notions like d-dimensional space-time. In recent 
years ambitious projects appeared where techniques are used to calculate physical processes in 
highly automatized ways HI SO. 

Here we focus on the calculation of Feynman integrals in Euclidean space-time. Though 
they are, in general, not the ultimate physical objects, their knowledge is very useful. First, if we 
know them analytically, analytic continuation gives a way to transform them to the Minkowskian 
region, if needed so. Moreover, if we solve, somehow, Feynman integrals involved in a given 
physical process analytically, we can use the knowledge in the Euclidean region to check these 
solutions numerically. This has been done in numerous cases at the 2-loop level (e.g. massive 
Bhabha scattering HEEIIZIEIEI, QCD calculations OS Oil [13), but also at higher loop levels 
or in more general context lfT3l[T4l[T5l[T6TD . For structural studies of quantum field theory, the 
numerical calculation of Feynman integrals in Euclidean space-time has also been proven to be 
useful, e.g. for the check of some conjectures in super- Yang-Mills theories fl 1 7L 1 1 811 . 

In the present article we describe two publicly available computational tools based on M at he- 
rn at ica which facilitate numerical calculations of the kind described. 

The first tool is the extended version v.2.0 of the AMBRE program [19] which generates Mellin- 
Barnes (MB) representations for Feynman integrals. We discuss the construction of MB-integrals 
with numerators of arbitrary rank R for L-loop cases. We also shortly report on additional features 
of older versions (v. 1.1 and v. 1.2) which were released after publication of H20L We explicitely 
work out a variety of non-trivial numerical examples. For this purpose, AMBRE is being combined 
with the Mathematica packages MB 11211 and MBresolve 112211 . 

The second tool is the Mathematica interface CSectors to the Ginac package sector_decom- 
position H23L The package sector_decomposition uses the sector decomposition method to 
calculate general polynomial structures which are present in calculations of Feynman integrals. 
In the spirit of AMBRE, we perform in CSectors for a given Feynman integral the automatic 
calculation of the characteristic F and U polynomials, add some normalizing factors consist- 
ing of Gamma functions and, finally, use a general formula for multi-loop tensorial Feynman 
parameterisations. The result is a user-friendly interface to sector_decomposition for the 
specific purpose of tensor Feynman integral calculations. What remains to be done by the user 
of CSectors/sector_decomposition is writing in a proper way the definitions of propagators 
(plus numerators, if they are present). Further, optional algorithmic strategies have to be chosen 
which are part of the core program [12311 . This is the stage of automatization reached also with 
AMBRE. 

The article is organized as follows. In Section [2] we prepare expressions for the general 
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multi-loop Feynman integral. Their evaluation based on Mellin-Barnes representations with 
AMBRE is described in Section [3j In Section [4] some details on sector decomposition and of using 
CSectors are discussed, and Section [5] contains numerical examples and few comparisons of the 
two approaches. It follows the Summary. In the appendix we list the most important Mathematica 
functions of AMBRE and options for CSectors. 



2 Definitions 



Comprehensive overviews of the presentations of Feynman integrals may be found e.g. in H24U25L 
Here, we repeat some basic formulae in order to define our notations. The L-loop Feynman integral 
in d = 4 — 2e dimensions with N internal lines with momenta qi and masses m ; , and E external 
legs with momenta p e is defined here as follows: 



G L [T R (k)} 



1 



d d k x ...d d k L T R (k) 



(in d / 2 ) L J {q\ -ml)*... (q 2 - m 2 )^ ...(q 2 N - m 2 )^ ' 
The numerator T R (k) is a tensor of rank R in the integration variables: 



and 



T R (k) 



Pi 



2 2 

q i - m t 



E 
e=\ 



Y^ilk-Pi 
1=1 



We allow for arbitrary indices v,-, the powers of propagator functions 

Next, the momentum integrals are replaced by Feynman parameter integrals: 



G L [T R (k)] 



-1 



r(vi)...r(%)i t=\ 



N 



N 



Y[dxixJ l 5(1 ~Y x d 



;=1 



jjN v -i(L+l)-R 
F N v -{h 



X 



(2.1) 
(2.2) 

(2.3) 
(2.4) 



f-4rrr(^V v -JL/2-r/2) F h[^.,^r \ {%S) 
where A^ v = I^i V;. 

The two functions U and F are characteristics of the topology of the Feynman integral. One 
may derive them from 



N 

JV = J>A = k(M L )k-2kQ + J, 
(=l 



(2.6) 
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where 

(M L )u> = J^auaa*i, (2.7) 



N 



i=l 
N 



Qi = J^aaPiXi, (2.8) 

N 

J = £(/f-m?)xf, (2.9) 



i=\ 

namely: 

U L {x) = det(M L ), (2.10) 
F l (jc) = -dot(M L )J + QM L Q, (2.11) 

with 

M L = det(M L ) {M L )-\ (2.12) 
The object s^ r ^R- r contains the tensor structure due to its two elements: 

*4> = 1, (2-13) 
^ = foroddr>0, (2.14) 

Sif-V-V* = gMlM2...|Mr-lMr 5 (2.15) 

and 



= 1, (2.16) 

^Ml-Mr = ^Ml...^Mr 5 (2.17) 

where we left out in the notations the indices related to the loop numbering, because they are fixed 
by ( |2.2J ) when the Lorentz indices are defined: 

~nm = (m l ~ 1 ) gnu, (2.18) 



L 

9>» = £ (M L ) U Q*. (2.19) 
1=1 



The product f^g/}^ 1 ' "'^ 1 ' '^j^'' "'^ j is completely symmetrized in its Lorentz indices; take as 
an example s&i^l'- 

+ st™ ^ 2 vA + j*£ p 9>f- + ^ p &>g v , (2.20) 
and, more explicitely, to e.g. T{1^ l k^) correspond the terms: 

sfa&^to = pVipVi^ (2.21) 
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and, with a different numerical factor (see (2.5 )): 



(2.22) 



For one-loop integrals, the rotation matrix M\ in the loop mom enta becomes trivial, M\ 
Mi = 1, and so also U\ = det(Mi) = 1, and Fi(x) = -J + Q 2 ^\The ( |2~5| ) then becomes: 



Gi(T R ) 



In case of e.g. L 



T(k^k^)^T(N v -d/2) 



f 1\N V r N N R W d+r\ 

11/=1 A l V *J J !=1 y=l r=0(-2)2 J F" 2 

(2.23) 



1,7? 



2, we get for the sum in (2.5 ): 



1 



2{N v -d/2-\) 



Fg 



M1M2 



(2.24) 



The general expressions as well as the examples agree with H19II . For one-loop tensors, equivalent 
expressions are also given in H26L 

An important observation is that the Feynman parameter representations for arbitrary Feynman 
integrals would be just dependent on polynomials in the x ? -, i.e. be sums of monomials in the %i 
with integer exponents, if there were not the two types of terms U(x) A ^ L,N,d,R ^ and F(x) B ( L,N,d > r \ 
The functions U (x) and F(x) are such polynomials, but they have non-integer exponents. While, 
the additional terms arising from the tensorial structure of the L-loop Feynman integrals are 
polynomials in the x,j^]The function U (x) = Y,n m n(xi) depends only on monomials m n (xj) and 
is positive semi-definite, while function F(x) — Y,n'[~ s n') m n'( x i) + U ( x ) x j m j depends also 
on the kinematical invariants and on the masses of the problem. For Euclidean kinematics, all 
the [—s n /] > and also F(x) becomes positive semi-definite H25II . One typical example is the 



F-function (3.4). 



The above formulae may be used for automated evaluations of specific Feynman integrals. For 
that purpose, one has to develop methods for their proper treatment, and two of them are worked 
out here. 



3 Integrations using Mellin-Barnes representations 

Iterated applications of Mellin-Barnes' formula 

1 1 



[A(s n )m n (x) +B(s n r)m n r(x)] a 2ni 



do[A(s n )m n (x)} a [B(s n i)m n/ (x)} 



r(fl) 



1 Sometimes it is useful to rewrite F\ (x) — > — + Q 2 , which agrees under the integral in ([23} with F\(x), 

but is now a bi-linear function of the X{. 

2 Because eqn. (13) of [ 19 1 was not sufficiently simplified, this was not evident there for L > 1. 
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may be used to transform the x- integrand of ( |2.5[ ) into a sequence of monomials in the Xj, allowing 
thus to perform the ^-integrations applying 

[ g^r - igf 

One remains with the problem to evaluate multi-dimensional complex Mellin-Barnes integrals. 

The publicly available package AMBRE prepares these Mellin-Barnes integrals. Since the first 
publication of AMBRE in H19I . several features has been added to the package. In AMBREv.1.1, the 
MB -representations can be constructed only for tensor integrals where all momenta of the tensor 
T[ki) are multiplied by external momenta, i.e. are part of scalar products. Since AMBREv.1.2 it is 
foreseen to generate MB -representations for just tensor integrals. Additionally, some new options 
were added, consult for details on them the webpages [20], where also appropriate examples are 
documented. One of the options allows to generate Feynman parameter representations without 
performing the ^-integrations, leaving them to be performed by the preferred technique of the 
user. 

Here we focus on AMBRE v.2.0 which generates tensor MB-integrals in a fully automatic 
way. Previous versions have the option Fauto which allows for manual manipulations on the F 
polynomial. Sometimes this is useful and helps to obtain a smaller dimensionality of MB-integrals 
(see e.g. IIT91 and the discussion for pentagons there). However, for many real processes with a 
large number of amplitudes, complete automatization is necessary. A second goal of the present 
version is a construction of MB -representations for tensorial planar n-\oop tensors. Though 
non-planar topologies can be also obtained directly with the loop-by-loop method employed in 
AMBRE 112711 . results may come out wrong j^] We recommend not to use AMBRE without independent 
checks for non-planar problems with several kinematical scales; they hopefully shall be treated 
more properly in the future. 

In order to include tensor structures of Feynman integrals, we had to modify properly the 
iterative procedure of H19H . The best way to explain this might be an explicit example, see Figure 

o 

Here all propagators are assumed to be massless and the numerator is of rank R = 3: 
{h. p ){k,.p){k 2 .p) 



[(* 2 - h ) 2 ]" 2 [(*! + P) 2 ]" 3 [k\] n * [(k 2 + p) 2 h 



The calculation starts by working out the sub-loop integration over k\ , which leads to the 
following F polynomial: 

F = [-k\] xix 2 + [-s] XiX 3 + [-(k 2 +p) 2 } x 2 x 3 . (3.4) 



According to (2.23 ), the generated tensor structure for the first sub-loop is an expression with 



five separate parts; we omit here a normalizing factor and end up with a table of five different 



3 We applied AMBRE successfully to e.g. the massless 2-loop vertex first studied in 11281 . which is a one-scale 
problem. A typical multi-scale problem is discussed in Section 3.2 of | 
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Figure 3.1: A simple example for the discussion of the tensor algorithm of AMBRE. 



integrals: 



-> (k^X 2 -p^X 3 )(k^X 2 -p^X 3 : 

-> {k^k^xi -kgpHxm, k^p^x 2 x 3 ,p^p^xl g^ 2 }. (3.5) 

Evidently, we will have to perform the second sub-loop integration over k 2 separately for all the 
above parts because they have different tensor ranks. As a rule, the rank of a given integral in 
the next step will include higher rank tensors than the original one. The situation after the first 
momentum integration can be symbolized as follows: 



/ 



X 



d d k 2 pimPiimfa-p) 



[kl]-^ [(k 2 + p)2]-3+e+ni+n 2 +n 3 + zl +Z2 



{k^k^ 2 x MBi, -k^p^ x MB 2 , k^p^ 2 x MB 3 , p^p^ 2 x MB 4 , g^ 2 x MB 5 }, 

(3.6) 



where the expressions MB, stand for different Mellin-Barnes parts of the net integral. These 



parts are different due to different Feynman parameters in ( |3.5[ ). As one explicite example we 
reproduce MB 1 : 

MBi = (-l) 2 - £ - Z2 (-s) Z2 r(2-£-n l -n 2 -zi)r(-zi)r(4-£-n l -n 3 -z2) 
x r(-z 2 )r(ni +zi + z 2 )T(-2 + £ + n l +n 2 + n 3 +zi +z 2 )/ 
[r(ni)r(n 2 )r(6 - 2e - m - n 2 - n 3 )T(n 3 )} . (3.7) 

Working with more than two loops, additional iterations can produce further 'fragmentations' of 
the expression. 

In the program, all above steps are hidden to the user. The only action is to define the input 
object: 
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invariants = {pl~2->s}; 

MBrepr [{kl*pl , kl*pl , k2*pl} , {PR [kl , , nl] *PR [k2 , , n2] * 

PR [k2-kl , , n3] *PR [kl+pl , , n4] *PR [k2+pl , , n5] } , {kl , k2}] 

The output is: 

{- ( ( (-1) " (nl+n2+n3+n4+n5) * (-s) " (4-2*eps-nl-n2-n3-n4-n5) *s~3 
Gamma [2-eps-nl-n3-zl] *Gamma [-zl] *Gamma [5-eps-n2+zl] 

*Gamma [4-eps-nl-n4-z2] *Gamma [4-2*eps-nl-n3-n4-n5-zl-z2] 

*Gamma [-z2] *Gamma [-4+2*eps+nl+n2+n3+n4+n5+z2] *Gamma [nl+zl+z2] 

*Gamma [-2+eps+nl+n3+n4+zl+z2] ) / (Gamma [nl] *Gamma [n3] 

*Gamma [6-2*eps-nl-n3-n4] *Gamma [n4] *Gamma [n2-zl] 

*Gamma[9-3*eps-nl-n2-n3-n4-n5-z2] 

*Gamma [-2+eps+nl+n3+n4+n5+zl+z2] )),...} 

(3.8) 

The dots stand for the remaining four MB -integrals. The complete expression can be found in the 
file MB_SE510m.nb at the webpage Il20ll . 



4 Integrations by sector decomposition 

The second approach - performing sector decompositions - transforms the x-integrand into a 
sequence of expressions where the singularities at d = 4, as a function of £ (d = 4 — 2e), are 
separated such that the arising expressions can be smoothly integrated. An appropriate algorithm 
for the automated computation of the £ series of multi-loop integrals has been formulated in 
lf3~0ll3Tll . In (231, a publicly available program is described which calculates the Laurent expansion 
of the following type of parametric integrals in £: 

d»x 8(1 - x» (f\xr eb ) ri [pM cj+edj ■ (4-D 

,j> i=1 V' =1 n=i 

The cii,bi, cj and dj are integers, and the Pj are polynomials in the variables x\, x n . The 
program may handle a product of several polynomials, and it is not required that the polynomials 
are homogeneous. These features are important for applying ( |4. 1) to our tensorial structures: it 



allows to calculate Gl(T) ( |2.5| ), and some examples are given in [23] for scalar integrals. 



Here we present an interface which directly calculates (2.5 ) using as the backbone the already 



programmed structure (4.1). What is expected by interface from the user is a proper definition of 
the integral ( |2.1[ ) and of the kinematical invariants to be used. Formally, it has to be done in the 
following way: 

invariants = {invariants as a rules}-; 

DoSector [{numerator}, {denominator}, {internal momenta}] [low, des] ; 
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Figure 4.1: Topology of the three-loop vertex exemplified in the text. 



(4.2) 



The first three entries in "DoSector" and "invariants" in (14.31) have the same form as in 



AMBRE. The "low" is the leading singularity of the considered integral. Sometimes it is not easy 
to determine this power of singularity. In order to be save, we propose to start with one or two 
powers of £ less than being relevant, giving zeroes for these cases. The "des" is the highest 
power of £ to which we are going to calculate an integral. One has to be careful because it is not 
possible to calculate only one £ term if it is not the leading singularity: If we want to calculate 
an e" term, but if lower terms in £ are present, we must start from the leading term. This feature 
arises because the pro gram calculates the whole Laurent series in £, taking into account also the 



pre factor jj^j^ m (2.5). 



The user input for CSectors is extremely simple. Here we reproduce a (slightly shortened) 
massless 3-loop vertex example with rank R = 3; the topology is shown in Figure [4~T 



Input: 



(4.3) 



<< CSectors. m; 
x=-ll; 

invariants={pl*p2->l/2*x,pl*p3->l/2*x,p2*p3->l/2*x,pl~2->x, 

p2~2->0,p3~2->0}; 
DoSectors [{k2*p2 , k3*pl , k3*pl} , {PR [k2+k3+p2 ,0,1] *PR [kl-k3 ,0,1]* 

PR [k3 , , 1] *PR [kl+k2-pl+p2 ,0,1] *PR [kl , , 1] *PR [k2 , , 3] } , 

{kl,k2,k3>] [-5,0] ; 
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Output: 



(4.4) 



CSectors by K.Kajda and V.Yundin ver:1.0 
last modified 22 sep 2009 

Using strategy C 
U & F polynomials: 

U = x4 (x5 x6 + x3 (x5 + x6)) + xl (x4 x5 + x2 

(x3 + x4 + x5) + x4 x6 + x5 x6 + x3 (x5 + x6)) + x2 
((x4 + x5) x6 + x3 (x4 + x5 + x6)) 

F = 11 (x3 (x2 + x4) + xl (x3 + x4)) x5 x6 

Qll = ... 
Q12 = ... 

Q21 = 121*x2*x3~2*x6~2/4 
Q22 = ... 

Generating C++ source . . . Intll . . . Intl2 . . . Int21 . . . Int22 . . . done 
Compiling source code ... Intll ... Intl2 ... Int21 ... Int22 ... done 
Running binary file Intll ... Intl2 ... Int21 ... Int22 ... done 

InputForm= 

{-7.622574999999999 - . 0260407/eps~4 
+ 0.049527000000000015/eps~3 - 0.4168788/eps~2 + 0.56955/eps, 
{1 . 3667737639753552 , 2 . 85804370185272*~-6/eps~4 , 

. 00009220935574625821/eps~3 , . 0004811810295896961/eps~2 , 

. 006549529654501916/eps}> 



The functions F and U in (4.4 ) are calculated with Mathematica, in the same way as in AMBRE. 



The parameters Qll ... Q22 are polynomials in the xi which arise due to the numerators in (2.5 ) 
for tensor integrals. These terms are select as positive semi-definite pieces Qll ... Q22 in order to 



satisfy the conditions of the sector decomposition algorithm. The dots for Qll, Q12, Q22 in (4.4 ) 
indicate expressions of some length and their explicit form can be found in file output_mercedes 
at the webpages lf32l . 



After "Generating C++ source" and "Compiling source code" in (4.4), the C++ files 



are running. They have the structure discussed in [23J and a user can inspect them by switching 



in the option "TempFileDelete->False" in (4.2). Other useful options are listed and described 
in the Appendix. 
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By default, the numerical errors of the results are also given. They are calculated by taking 
the combined error for all the integrals /, calculated at a given term of the e-expansion: 



AI F 



\ 



(4.5) 



As it is well known, some of integrals, especially massless ones, can be difficult to integrate 
numerically. Here the proper choice of a sector decomposition strategy may help H23L In ( 4.4[ ), 
the calculation has been done with optional strategy C. In order to control which e-term is currently 
calculated, the default option TempFileDelete->False produces a log file. For the Qll term of 
the example, it is: 



» SE510m_numll.log « 
Order eps~(-2) : +/- 
Order eps~(-l) : +/- 
Order 1: -203.056 +/- 0.00948256 
Order eps: 153.241 +/- 0.0721185 



A list of all options can be obtained with the command ?0ptions in Mathematica. 

In order to test CSectors we have checked many topologies: multi-loop tadpoles and self- 
energies, three-loop vertices up to rank R = 5 and with double dots on propagators (corresponding 
to setting index v = 3, because a 'dot' raises the index by one), four-point functions up to 
rank R = 4, and some one-loop five- and six-point functions. For some higher rank tensors we 
have used Integration-by-Parts decompositions of the Feynman integrals using the computer 
algebra package IdSolver^J Some numerical examples of these tests can be found in the file 
numerical_checks .nb at the webpages 1321 . 



5 Numerical results 

In Section[3]it has been shown how to define Feynman integrals and how to get numerical values for 
them at chosen kinematical points using the MB-method. For sector decomposition, the same has 
been discussed in Section]?} For scalar integrals it is straightforward to use, together with AMBRE, 
the MB package, and to perform the numerical integrations [12 1 11 . For tensor integrals, especially 
with loop order L > 1, we have usually many MB-integrals for which the command MBrules of MB 
has to be performed in order to find proper integration paths for the MB -integrations. Sometimes, 
depending on the degree of divergency of the original Feynman integral, this is not possible if 
only an analytical continuation in e is done. Then, an analytical continuation in one of the indices, 
may be successful. This has also been automatized. If MBrules does not find a valid rule for some 
£ ^ 0, then the power of the first propagator Vi is changed, e.g. if Vi = 1, then Vi — > Vi = 1 + r\ 

4 IdSolver is an unpublished C++ package performing Feynman integral reductions with the Laporta algorithm. 
J.G. thanks M. Czakon for the opportunity to use it. 
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is applied. If again MBrules cannot find a valid rule, V 2 is treated analogously, and the procedure 
can be continued until MBrules is succes sfulS 

We introduced the auxiliary file MBmim.m to MB.m which realizes this procedure and the 
subsequent automatized analytic continuation, e-expansion, and numerics. The file may be 
obtained from the webpages H20L 

If an already prepared MB-representation repr is available, e.g. in (3.8 ), it is enough to use 
the MBnum function (see Appendix for details), e.g.: 

MBanalytic=MBnum[repr , -1, {s -> -11}, 

{nl -> 1, n2 -> 1, n3 -> 1, n4 -> 1, n5 -> 1}, 2] 
res=MBintegrate [MBanalytic , {s->-ll}] 

This gives the following numerical result (see also file out_SE510m at the webpages f20l0 : 

SEnumMB = -7.5625/e 2 - 20.4506/e - 178.18 ±0.0171936 + (18.3642 ±0.0248465) e 

(5.1) 

The corresponding values returned by CSectors are (see also file output_SE510m at H32II ): 
SEnumSD= 

{-178.1927 - 7.56258/eps~2 - 20.4505/eps + 18 . 394000000000005*eps , 
{0 .011130644628528934, . 00029563/eps~2 , 
. 00302109/eps , . 06629380324170578*eps} 
> 

(5.2) 

Another example is two-loop Bhabha scattering. For two-loop 4-point functions, sector 
decomposition needs a lot of RAM (typically up to few GB) and also of computing time. In such 
cases, often the numerical integrations are done faster using the MB-method. Some numerical 
results for the Bhabha Feynman integral with rank R = 2, 

Bl = [d%d% ik l .p l ){k 2 . P2 ) 



{k\-m 2 ){h + Pi) 2 [(h +P1 + P2) 2 ~ m 2 } 

1 (5.3) 



(h - k 2 ) 2 (k 2 - m 2 ) [(k 2 +pi+ p 2 ) 2 - m 2 ] (k 2 + pi+p 2 + p A ) 2 ' 

corresponding to the topology shown in Figure 5. 1[ can be found in Table |5 . 1] 

The Mandelstam variables are s = (p\ + p 2 ) 2 and t = (p 2 + P4) 2 . The so-called reducible 
numerators of a tensor Feynman integral can be contracted with propagators. Apart from speeding 
up calculations, such contractions have been used to check the implementation of tensor structures 
into AMBRE and CSectors. In our example, the relation 

(h-pi) = \[{k l +p l )- 2 -{k 2 -m 2 )- l -2m 2 ] (5.4) 



5 If it happens that an analytical continuation in only one additional parameter 77 is not sufficient, the program will 
stop with a proper remark. 
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Figure 5.1: The planar two-loop Bhabha topology Bl. 
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m = 1 




0.65734(6) 


0.659(2) 


3.186(6) 


3.174(2) 




0.13921(8) 


0.1396(5) 


1.0383(1) 


1.0381(3) 


£ 2 


0.018095 


0.01835(9) 


0.28817(1) 


0.28816(4) 


£ 3 


0.104974 


0.10500(2) 






£- 4 


-0.0217857 


-0.021785(3) 






T M 


368 


26700 (7.5 h) 


945 


70220 (19.5/i) 






s = —5, t = 


= -7 





Table 5.1: Numerical values for the first terms of the £ expansion of the massless and massive 
Bhabha 2-loop double box with tensor rank R = 2, defined in ( |5.3[ ). The topology is shown in 
figure |5TTj The package CSectors was used with strategy X for m = and with strategy C for 
rn^O. 



may be used to change the rank R = 2 tensor of ( |5.3[ ) into three tensor integrals of rank R = 1 (in 
addition reducing the number of propagators by one for two of them). In this way, the numerical 
results given in Table |5. l] can be cross-checked. 



Further, let us present a four-loop self-energy scalar diagram, Fig. |5.2[ With CSectors it 
takes a few hours to calculate the constant term of the £-expansionjj As it is a scalar diagram, 
we can make direct comparisons with FIESTA2 03H . which is much faster but applies to scalar 
integrals only. For AMBRE/MB the complete calculation takes about two minutes only, which is 
comparable to using FIESTA2. 

Finally, as a more complicated case, let us take a pentabox of rank R = 3: 

{h-h){h-k 2 ) (h-k 5 ) 



PB 



d d hd d l 2 



1 



X 



(h - h - k 2 - k 3 ) 2 ll(l 2 - k 5 f(l 2 - k 4 - k 5 ) 2 (h + hf 



(5.5) 



6 More precisely, it takes eight ours both with strategies C and X. This and the other examples calculated in this 
paper have been run on a Xeon personal computer. 
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Figure 5.2: Four-loop self-energy. For details of its calculation see the files MB_SE41oop.m, 
MB_SE41oop . out and SD_SE41oop . sh, SD_SE41oop . out at the webpages 112011321 . 



Figure 5.3: The pentabox topology PB. 

Its topology (called PB here) is shown in Figure [573] 

All external momenta k{ are incoming, sij = (ki + kj) 2 , S12 = — 1,^23 = -2,534 = —3,545 = 
—4,5i5 = — 5), and the numerical result is: 

PB= 

{-7 . 431348973217098+0 . 328125/eps~4+0 . 5340786885498234/eps~3 

-0.44570891627426246/eps~2 - 3 . 1058689125651284/eps , 
{0 .012739775270198513+5 . 560657237081255*10' (-10) /eps~2 

+0 . 0003572744410650107/eps}} 

(5.6) 

Some caution must be paid to the higher dimensional MB integrals. They can give underesti- 
mated errors, taken from the CUBA error output. Whether this happens can be checked switching 



on the MB.m option Debug. The result in ( |5.6| ) has been obtained for the following set of MB 
parameters (see the file MB_PBox .m at |[20l ): 

SetOptions [MBintegrate , Verbose->False, PrecisionGoal->10 , 
AccuracyGoal->24, MaxPoints -> 1000000] ; 



15 



As it can be seen in MB_PBox . m, the construction of the MB -representation for ( |5.5[ ) starts with 
momentum I2, leading to maximally nine dimensional MB -integrals. Applying Barnes lemmas 
reduces the integrals to at most seven-dimensional. Starting with the integration over l\, we get 
maximally eleven-dimensional integrals. Altogether, the calculation takes about five hours: The 
generation of the MB-representations needs a couple of minutes, and the rest of the time is needed 
for analytic continuations and numerical integrations. 

Additional instructive examples for using AMBRE and CSectors may be found at the webpages 

mm. 



6 Summary 

The CSectors package has been prepared as an interface to the sector decomposition package 
sector_decomposition for the automatic numerical evaluation of tensor Feynman integrals. 
The AMBRE package has been extended for the automatic treatment of tensor structures of multi- 
loop Feynman integrals. 

For CSectors, the bottleneck is the numerical evaluation based on sector_decomposition 
and Ginac which, especially for higher rank tensors and higher loop orders, consumes a huge 
amount of RAM. However, for smaller problems, up to two loops and tensors of moderate rank, 
the program works well. For more complicated cases we recommend to use the Mellin-Barnes 
approach. Here there are many possibilities to optimize the way to get numerical results. This 
can be done at the level of construction of MB-representations (e.g. by a change of order of 
the integration over internal momenta), then there are different ways of analytic continuation, 
and finally the Mathematica package barnesroutines [13411 can be used to try to reduce the 
dimensionality of the MB-integrals; examples for the latter can be found in H35H and at H20[ 13210 . 

For the near future it is planned to automatize the construction of MB-representations for 
non-planar Feynman integrals, what deserves to leave the loop-by-loop approach. Further, it is 
foreseen to build MB-representations for special forms of linear propagators which appear e.g. in 
heavy quark effective theory (HQET) or in calculations of the QCD static potential. 
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7 Appendix 

The instructive examples for using AMBRE and CSectors may be found at the webpages H20M32L 



AMBRE (ver l.X) 

The appropriate loop integral is defined by: 

(1) a list of kinematic invariants, e.g. invariants = {pl*pl -> s}, 

(2) Fullintegral[numerator, propagators, internal momenta, options], 
invariants must be defined before Fullintegral. 

The arguments of Fullintegral are as follows: 

• numerator : numerator which can be given in the contracted form, e.g. {kl*pl ,k2*p2> 
or in the uncontracted form, e.g. {k2 [mul] ,k2 [mu2] } , 

• propagators: product of propagators of the form PR [q,m,nl] = (q 2 — m 2 ) ni , 

• internal momenta: list of internal momenta, e.g. {kl ,k2}. 

This version of AMBRE uses a semi-automatic approach when building Mellin-Barnes repre- 
sentation. That methodology is accomplished by the following two functions. 

(3) IntPart[iteration, options] - prepares a subintegral for a given internal momentum by 
collecting the related numerator, propagators and integration momentum: 

• iteration: iteration for which subintegral will be prepared. In practice IntPart function 
must be executed in specific order i.e firstly IntPart [1] then IntPart [2] and so. 

• options: 

- Text: it can have two boolean values True or False. Controls if additional text 
appears during calculations. 

(4) SubLoop [integral, options] - determines for the selected subintegral the U and F polynomials 
and an M-B representation. 

• integral: this argument must be left as it is. 

• options: 

- Text : as in the IntPart function. 

- Xintegration: controls whether integration over Feynman parameters is performed 
or not. 
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(l)-(4) constituates basic functions. An additional functions are: 

(5) Fauto[value] - allows user specified modifications of the F polynomial f upc. Must be used 
after IntPart and before SubLoop. 

• value : can be or 1. For the first one user can modify F polynomial. For the latter this 
possibility is turned off. 

(6) BarnesLemma[representation,number,options] - function tries to apply Barnes's first or 
second lemma to a given representation. 

• representation: M-B representation to be checked. 

• number: it has two possible values, 1 for the first Barnes lemma and 2 for the second 
lemma. 

• options: there are the following two boolean options for this function 

- Text : displays or does not display an additional information. 

- Shifts : searches for the pairs of two integration variables Zi + Zj and Zi — Zj which, 
after application of the appropriate shift one of it is cancelled. 

(7) ARint[result,i_] - in version 1 . 2, it displays the MB -representation number i for Feynman 
integrals with numerators 

AMBREnLOOP (ver 2.0) 

The basic functions of AMBREnLOOP are: 

(1) MBrepr[numerator, propagators, internal momenta, options]. 

It returns M-B representation for a given loop diagram. All arguments, except options, have the 
same form as parameters in version 1 . X. 

• options: 

- Text : displays or not information text. 

- OptimizedResult : the final MB representation is written in such a way that Gamma 
functions are factorized. This option makes an output more condense. However sum 
of different Gamma functions for a given integral can cause problems if treated as it 
stands in analytic continuation (finding rules). Optionally it is switched off. 

- BarnesLemmal : turns on or off first Barnes lemma checking. 

- BarnesLemma2 : does the same as above option but for second Barnes lemma 

Intermediate output of a given subloop is displayed in a specific format using INT function: 
INT [numerator , representation, propagators 1 ,propagators2] 
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• argument s : 

- numerator : it is, for tensor integrals, of the form {k2 [mul] , k2 [mu2] }. 

- representation: Mellin-Barnes representation obtained during the calculation of 
the current sub-loop/iteration part. It is multiplied by Mellin-Barnes representations 
which were calculated in the previous step/iteration. 

- propagators 1 : propagators extracted out of the F -polynomial in the current itera- 
tion. 

- propagators2 : keeps propagators in the same form as propagatorsl does. The 
propagators2 are independent of the present integration variable but undergo next 
iteration(s). 

(2) BarnesLemma[representation,number,powers of propagators,options] - this function 
works as in the previous versions of AMBRE. The only difference is the extra parameter: 

• powers of propagators: - here a list of powers of propagators which appear in an input 
loop integral must be given, e.g. {nl , n2 , n3>. 

CSectors 

The appropriate loop integral is defined by: 

(1) a list of kinematic invariants, e.g. invariants = {pl*pl -> s}, 

(2) DoSectors[numerator, propagators, internal momenta, options] [min, max] 

The arguments of this function are exactly the same as in case of MBrepr in the AMBREnLOOP 
package. However, in contrast to the AMBRE package, numerator can be written only in the 
uncontracted form. 

The package allows to modify its behaviour by adding additional options: 

• SetStrategy: chooses one of the strategies available in sector decomposition libraries 

E3. 

• SourceName : a prefixing for source, binary and log files; the option just allows to choose 
any name for the files connected with the calculation of a given integral. 

• TempFileDelete : by default it is set to TempFileDelete->True; when set to False, it 
does not delete C++ source and binary files as well as log file. 

• LogFile : causes CSectors to create (or not) a log file, where numerical results for given 
integral and epsilon are stored. The default is True 
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• ShowErrors : controls whether the errors of the numerical calculation will be displayed or 
not; errors are calculated using the function res . get_error () of sector decomposition 

E3. 

• IterationsLow, IterationsHigh, CallsLow, CallsHigh: these are Monte-Carlo 
parameters, see lT2~3Tl for a description. 

• compiler: allows to choose another compiler; the default is g++. 

• cppflags, libs : the paths to header and library files, required by sector decomposition 
libraries IIZ51 . 

• min , max : indicates minimum and maximum of the Laurent series expansion in epsilon. 
The default options can also be displayed by the command Options [DoSectors] . 
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